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ABSTRACT 

An attempt is made to determine the three-dimensional ocean circulation from satellite altimeter measurements 
by assimilating Geosat sea surface height data into an eddy-resolving quasigeostrophic (QG) model of the eastern 
North Atlantic Ocean. Results are tested against independent information from hydrographic field observations 
and moored current meter data collected during the Geosat ERM. The comparison supports the concept of 
inferring aspects of the three-dimensional flow field from sea surface height observations by combining altimetric 
measurements with the dynamics of ocean circulation models. 

A Holland-type QG model with open boundaries is set up on a 2000 km X 2000 km domain of the eastern 
North Atlantic between 25 and 45°N, 32° and 8°W. By using a simple nudging technique, about two years of 
Geosat altimeter data are assimilated into the model every five days as space-time objective analyses on the 
mode! grid. The error information resulting from the analysis is used during the assimilation procedure to account 
for data uncertainties, Results show an intense eddy field, which in the surface layer interacts with a meandering 
Azores Front, Compared to Geosat, the model leads to smoothed fields that follow the observations. 

Model simulations are significantly correlated with hydrographic data from March 1988 and June 1989. both 
close to the surface and in the subsurface. Good agreement is also found between the model velocity fields and 
moored current meter data in the top two model layers. The agreement is visually weak in the bottom layer 
although a coherence analysis reveals an agreement between the model simulation and current meter data over 
the full water column at periods exceeding 80 days. 


1. Introduction 

During the last decade satellite altimetry has become 
a powerful tool for observing and understanding the 
ocean circulation by providing global and continuous 
sea surface height (SSH) observations. It was shown 
that the altimeter signals even from Seasat and Geosat 
were capable of adequately representing ocean surface 
dynamics (e.g., Willebrand et al. 1990) and could be 
used to monitor the surface current fields (Gordon and 
Haxby 1990; Stammer et al. 1991). 

However, it is not obvious how to reconstruct the flow 
field of the full water column, given only surface flow 
observations. The related problem of the observability 
of the 3D flow field from surface observations has lead 
to a controversial discussion and stimulated many in- 
vestigations. Observations show' that the vertical struc- 
ture of the ocean is predominantly low' mode (e.g., Miil- 
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ler and Siedler 1992; Fukumori and Wunsch 1991). 
which suggests that the interior ocean current field is 
effectively coupled to the sea surface. This is supported 
by results from numerical studies, which show that the 
abyssal flow field is constrained by the sea surface pres- 
sure boundary condition (e.g., Holland and Malanotte- 
Rizzoli 1989, hereafter HMR; Haines 1991 ). The com- 
bination of altimetric observations with numerical mod- 
els should therefore allow one to construct the flow field 
over the full water column consistent with surface ob- 
servations. 

To assimilate altimeter data into numerical circulation 
models, various assimilation techniques with a wide 
range of complexity are available. Ghil and Malanotte- 
Rizzoli (1991) provide a comprehensive review of as- 
similation methods presently under use in the atmo- 
spheric and oceanographic context, and Anderson and 
Willebrand (1989) review their oceanographic appli- 
cations. Inverse methods allow us to account for errors 
in both data and models and are all constrained and 
unconstrained model-data combinations. But for real- 
istic problems on large space or time scales, optimi- 
zation methods known to be useful, such as the Kalman 
filter smoother or the adjoint methods, overwhelm pres- 
ent-day computer resources. A few example applica- 
tions in the context of altimetric data can be found in 
Moore (1991), Schrdter et al. (1993), and Fukumori et 
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al. (1993). Instead, simple methods like nudging or di- 
rect insertion have been used widely in feasibility stud- 
ies, which aimed at reconstructing the interior ocean 
flow from surface data. Most of those studies are based 
on twin-experiments with quasigeostrophic (QG) mod- 
els in which synthetic observations are assimilated 
(Hurlburt 1986; DeMey and Robinson 1987; Berry and 
Marshall 1989; Verron and Holland 1989; Holland and 
Malanotte-Rizzoli 1989; Haines 1991). But extended 
nudging techniques have also been developed recently, 
to assimilate altimeter data into primitive equation (PE) 
models (e.g., Mellor and Ezer 1991 ; Haines et al. 1993). 
Similar procedures were applied to realistic observa- 
tions both with QG and PE models (e.g., Holland et al. 
1992; White et al. 1990; Robinson et al. 1989; Oschlies 
and Willebrand 1995). 

However, twin-experiments can be vastly unrealistic 
because synthetic observations are not contaminated by 
measurement errors and are perfectly consistent with 
the underlying model physics, both in its mean and eddy 
components. In addition, synthetic observations are 
available to arbitrarily high resolution in both space and 
time. These situations are much in contrast to realistic 
conditions. Altimeter data are provided only along sub- 
satellite tracks and are not fully compatible with the 
underlying model. In addition, they are plagued by mea- 
surement errors and data loss. Thus, they provide only 
a crude estimate of the true space-time variability of 
the ocean, which, in general, will differ from that of the 
model. Moreover, the model mean field will differ from 
the observed mean state. 

Under such circumstances, the ocean model has to 
play a much more crucial role in extracting the dynam- 
ical signal from the noisy data and in carrying the in- 
formation into data-sparse areas, both horizontally (be- 
tween subsatellite tracks and in data dropout areas) and 
vertically from the surface toward the interior ocean. 
Faced with these difficulties, it might be questioned if 
the earlier encouraging results from identical twin-ex- 
periments are applicable to the real world. Unfortu- 
nately, it was hardly possible to compare model results 
with in situ data, due to the lack of simultaneous sub- 
surface observations. 

This latter issue is addressed in this study, in which 
Geosat observations are assimilated into a regional 
eddy-resolving quasigeostrophic ocean circulation mod- 
el of the Iberian Basin using a simple nudging technique 
(Verron and Holland 1989). Geosat anomalies are ref- 
erenced to the mean dynamic surface topography pro- 
vided by Robinson et al. (1979, henceforth RBS). The 
model's success in simulating the true ocean circulation 
in subsurface layers is tested subsequently against in 
situ hydrographic and current meter data. 

In a parallel and much similar effort, Capotondi et 
al. (1995a,b) investigated the dynamical consequences 
of assimilating both time mean and time-varying com- 
ponents of sea surface height data in the Gulf Stream 
extension area. Mechanisms of the model adjustment 


and the characteristics of the resulting new mean state 
are investigated when both components are assimilated 
separately. It was shown that the model behavior in the 
presence of the surface constraint can be described in 
terms of vertical Foffonof modes. A prescribed mean 
surface component determines a distortion of the geo- 
strophic contours in subsurface layers, which constrain 
the deep mean flow field. The intensity of the mean flow 
is specified by the inflow conditions, as well as eddy 
forcing and dissipation. 

Many of those dynamical considerations regarding 
the constraining of the mean flow field by surface data 
are applicable to our study and will not be repeated here. 
Instead, the primary focus of this paper is to test the 
success of reconstructing the time-varying flow field in 
the full water column. Because the boundary current 
area is dynamically distinct from the eastern North At- 
lantic and because of differences in various technical 
details, we consider our approach in the eastern North 
Atlantic to be distinct from that of studies in the western 
Atlantic. Conclusions drawn from both studies are much 
in agreement, however. 

This paper is organized as follows. Section 2 ad- 
dresses problems related to regional circulation models, 
the model initialization, and the assimilation procedure. 
In section 3 results of the assimilation experiments are 
presented, and section 4 deals with the test of the model 
simulations using the independent information of hy- 
drographic and current meter data. 

2. The model and assimilation procedure 

a . Model of the eastern North Atlantic 

Our experiments are based on Holland’s (1978) QG 
model of a closed basin, which was formulated for N 
arbitrary layers by Chow and Holland (1986). This QG 
model has been widely used during various eddy-re- 
solving experiments with both closed and open bound- 
aries. Holland (1986) gives a summary of many of those 
studies. 

For the present study, the QG model was set up in a 
2000 km X 2000 km model domain that covers the 
eastern North Atlantic, north of the Canary Islands and 
east of the Azores between 25° and 45°N, 8° and 32°W. 
The resolution is 10 km in the horizontal and three layers 
in the vertical. A time step of 5400 s was used. In the 
bottom layer, 20% of the realistic bottom topography, 
truncated above 2000-m depth, was included. Model 
parameters used during our experiments are listed in 
Table 1, most of which are based on earlier studies by 
Holland and Malanotte-Rizzoli (1989). 

Figure la shows the real topography in the model 
domain. While the continental shelf is a natural bound- 
ary at the eastern side, the bathymetry is dominated by 
the Mid-Atlantic Ridge (MAR) in the west. In between, 
the ocean floor can be considered approximately level 
at a mean depth of about 4500 m. Various field studies 
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Table L Model parameters. 


L y x L, = 2000 km X 2000 km 

Basin dimensions 

<5, = 10 km 

Grid length 

8t = 5400 s 
W, = 500 m 

Time step 

H, = 1000 m 
H, = 3000 m 

Layer thicknesses 

4> t = 35°N 

Reference latitude 

A = 20°W 

Reference longitude 

f t = 8.365 X 10 ’ 

Coriolis parameter at refer- 
ence latitude 

P = 1.875 X 10 " 

Rate of change of Coriolis 
parameter with latitude 

A 4 = 8 X 10 9 m 4 s 1 

Biharmonic friction coeffi- 
cient 

€ = 1 X 10 7 S 1 

Bottom friction coefficient 

g|/> = 9.63 X 10 4 m s 2 

Stratification parameters 
g(A p/p (( ) 

= 1.95 X 10 ■’ m s ! 

At interfaces between lay- 
ers 


have shown (e.g., Kase et al. 1986, 1989) that the hy- 
drography in this area is dominated by a mesoscale eddy 
field superimposed on the permanent frontal structure 
of the Azores Current (AC). Using the hydrographic 
data of Robinson et al. (1979) from a 1500-dbar ref- 
erence, the mean AC is observed in Fig. lb to flow 
eastward at approximately 5 cm s _l , near 43°N. The 
Portugal Current is also indicated near 40°N by an en- 
hanced current speed east of 20°W. 

Using an AMayer model, the vertical structure of the 
(N — \) th baroclinic mode leads to the appropriate dis- 
cretization of the vertical layers (Flierl 1978). Figure 2a 
shows a mean cr H profile calculated from the RBS data 


over the model domain together with its corresponding 
profile of the Brunt-Vaisala frequency 


n = fj*. m 

Related profiles of the first three eigenmodes (barotropic 
and first two baroclinic) are displayed in Fig. 2b, which 
follow as a solution of the eigenvalue problem 


1 

dz y/V 2 ) dz 


+ A Z F = 0, 


(2) 


subject to the boundary condition dFIdz = 0 at z=0, 
-H, assuming a mean ocean depth of H = —4500 m. 
Related Rossby radii of deformation r. = A, l/2 , i = 0, 

l , 2 of r u = 2330 km, r, = 25.7 km, and r 2 — 10.3 km, 
respectively, are in agreement with Emery et al. (1984). 
Based on the structure of the second baroclinic mode 
(zero crossings are at 350 m and 1580 m, respectively) 
the depths of the three layers have been chosen as 500 

m, 1000 m, and 3000 m, respectively. Together with 
layer-averaged densities this results in radii of the first 
and second baroclinic model modes of 26 and 13 km, 
respectively. 

The regional nature of the model and the need to 
simulate the Azores frontal structure requires some 
thoughts on boundary conditions. To simulate only a 
fraction of the real ocean requires lateral boundary con- 
ditions that represent the influence of the large-scale 
basin circulation on the local model domain. We have 
used boundary conditions as suggested by Charney et 
al. (1950), which prescribe the model streamfunction 



Fig. 1. (a) Realistic bottom topography in the model domain of the eastern North Atlantic. Marked are the positions of two moorings 
from which current meter data are used for model verification (compare Fig. 14). (b) Mean geostrophic surface circulation in 50 m relative 
to 1500 dbar as inferred from Robinson et al. (1979) climatologic atlas data. 
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Fig. 2. (a) Mean profiles of potential density anomaly (T h and Brunt-Vaisiila frequency N derived from the Robinson el al. (1979) hydro- 
graphic data for the model domain, (b) Barotropic and first two baroclinic vertical eigenfunctions. 


'P, and the relative vorticity £ k = V :x P A on the complete 
boundary in each layer k during inflow. During outflow 
conditions, £ is extrapolated from the interior model 
domain. Boundary conditions for 'P and £ are taken from 
the climatologic dynamic topography at the mean layer 
depths of the first and second layers (250 and 1000 dbar) 
determined from objective analyses of RBS atlas data 
relative to 1 500 dbar (Fig. 4a, b). [Details on the analysis 
of the climatology are provided by Stammer et al. 
(1991).] Due to the lack of detailed information about 
the bottom layer, it is assumed to be at rest initially with 
no water entering or leaving subsequently. 

In addition to lateral boundary conditions it is nec- 
essary to specify suitable surface boundary conditions, 
which also include the model forcing by the altimeter 
data. Generally, it is desirable to use a complete surface 
wind stress field taken simultaneously with the altimeter 
measurements. For example, daily ECMWF-model 
analyses would be a suitable source of wind data. But 
because of problems that arise in basically all models 
when simulating the Azores Front by real wind stress 
fields, we make the alternative assumption that the as- 
similation of altimeter data provides a vorticity source 
or sink at the surface equivalent to wind forcing, which 


is necessary to keep the model closely tied to the Geosat 
observations. 

b. Model initialization 

The relatively short periods over which real altimeter 
data exist and the regional nature of the model require 
a model initialization procedure that extrapolates the 
initial surface information into the deep layers in a dy- 
namically consistent way that also takes the climato- 
logical background field into account. Various such ap- 
proaches exist (e.g., see DeMey and Robinson 1987). 
Here we follow Holland et al. (1992), who used an 
initialization method based on the conservation of the 
QG potential vorticity U in the interior ocean. During 
this procedure, the initial climatological potential vor- 
ticity fields n* from layers k > 1 are inverted jointly 
with an initially observed field 'P l)hs in the surface layer 
to infer a dynamically consistent set of model stream- 
functions 'P,, 'P,, 'P* and of potential vorticity 11,, n>, 
II, in the entire water column. This initialization assures 
that the potential vorticity of water particles, received 
initially from the climatological background, is consis- 
tent with subsequent source of II during inflow con- 
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ditions. Haines (1991) uses a similar method to assim- 
ilate altimeter data into an ocean model through suc- 
cessive reinitialization. 

Figures 3a and 3b display the climatological stream- 
functions 'F,, of the mean ocean circulation for layers 1 
and 2, respectively, which have been inferred from RBS 
data based on the dynamic topography in 250 and 1000 
dbar relative to 1500 dbar. The fields are dominated by 
the large-scale current structures within the model do- 
main, with an inflow at the western boundary in both 
layers related to the eastward flowing Azores Current 
and to the extensions of the North Atlantic Current 
(NAC) in the northern part of the model domain. The 
resulting climatological potential vorticity fields FI,, for 
layers 1 through 3 are shown in Figs. 3c-3e, which are 
based on the assumption of a resting deep ocean. While 
the top layer is characterized by an enhanced horizontal 
FI gradient associated with the AC (Fig. 3c), the lower 
layers are dominated by the planetary potential vorticity 
(Fig. 3e). 

Figures 4a, 4b, and 4f display the initial model 
streamfunction fields 'FX?,,) in the three model layers, 
which result from the inversion of the observations ¥ obft 
(Fig. 4a) in the top layer and the subsurface climato- 
logical potential vorticity II,, (Figs. 3d and 3e). The 
resulting initial eddy field relative to the climatological 
background (Fig. 4d to 4f) shows a strong vertical co- 
herence with amplitude decreasing with depth by one 
order of magnitude and eddy scales increasing toward 
the bottom. During this type of model initialization, the 
surface eddy amplitude is projected predominantly onto 
the barotropic and first baroclinic mode, while the sec- 
ond baroclinic mode stays approximately unchanged. 

c. The assimilation procedure 

To assimilate altimeter data into the QG model, a 
simple Newtonian relaxation method (“nudging”) is 
used. This technique has been used in many previous 
altimeter data assimilation studies and is extensively 
reviewed by Verron and Holland (1989) and Holland 
and Malanotte-Rizzoli (1989). In order to nudge the 
model streamfunction toward altimeter observations, the 
vorticity equation of the surface layer is extended to 
include the relaxation term /?(V 2x F (lbs — V 2 ^,), where 
the QG streamfunction ' v F obs is related to the altimeter 
observations h ohs of the surface elevation as 

V*. = f A.*,- (3) 

The weight of the relaxation term relative to the model 
physics is determined by the timescale r r — R 1 and by 
the difference between the model and observations. 

Empirically determined coefficients from twin-ex- 
periments suggest that r r should be of the order of Vi 
to 1 day (see HMR). The latter value has been used in 
the assimilation experiments of Holland et al. (1992) 
during which a QG model was nudged toward Geosat 


observations in the Agulhas retroflection area. Capo- 
tondi et al. (1995a,b) used r t . = Vi day in the Gulf Stream 
region. A scale analysis shows that those numbers are 
about an order of magnitude smaller than all time scales 
of the remaining physical terms in the equation (see 
Capotondi et al. 1995a), which implies that the model 
surface streamfunction is effectively replaced by the ob- 
servations. 

However, real altimeter data are subject to a variety 
of errors that arise from such problems as incorrect en- 
vironmental corrections or the lack of an instantaneous 
ocean background field. To account for those errors it 
seems plausible to use a somewhat reduced weight on 
the data by diminishing the relaxation coefficient to al- 
low model physics to correct tendencies in the sea sur- 
face elevation that are not dynamically consistent. This 
issue will be addressed further below. 

During all experiments, the observations /r obs were 
provided as a space-time objective analysis on the mod- 
el grid every 5 days, estimated previously from the 
alongtrack data using correlation scales in space and 
time of A = 100 km and r = 15 days, respectively [see 
Stammer et al. (1991) for details on the objective anal- 
ysis of the Geosat data]. The absolute surface elevation 
/iohs — h + h r was obtained from the sum of the Geosat 
anomalies h f and the climatologic dynamic topography 
h from the center of the first model layer at 250-m depth 
relative to 1500 dbar. An example of the resulting fields 
is shown in Fig. 4a. 

To assimilate those altimetric fields, the relaxation 
coefficient R was modified to include a positive weight- 
ing function e, and a Gaussian exponential in time to 
be R = e R 0 exp[( — \tlrp) 2 ], where A/ = t - f obs is the 
distance between the model and the observations in time 
and t p is an exponential decay (persistence), which can 
be represented by the decorrelation time of the obser- 
vations. From Geosat observations we estimated r p to 
be of the order of 5 days. The coefficient e accounts for 
the Geosat estimation error provided by the objective 
analysis. It is close to unity where the uncertainty of 
the objective analysis is negligible, and it approaches 
zero when the analysis uncertainty is close to the process 
variance. 

Instead of using objective analyses, one could also 
directly assimilate the alongtrack Geosat data. This ap- 
proach has been successfully used by Verron (1992). 
Here we do not follow this method because it would 
constrain the model more by uncertainties than is the 
case when using the analyzed fields. This holds in par- 
ticular for the tidal alias error (Schlax and Chelton 1 994) 
by which an error signal is introduced with apparent 
westward phase propagation similar to first mode bar- 
oclinic Rossby waves. The tidal alias error is much re- 
duced in the space-time objective analyses used in this 
study. 

A schematic of the model setup is drawn in Fig. 5. 
The hatched area next to the boundaries indicates a fric- 
tion layer that is built in to damp reflected planetary 
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Fig. 3. (a) and (b) Mean fields of the climatological QG streamfum 
namic topography in 250 and 1000 dbar, respectively, relative to 
climatologic potential vorticitv IF, based on the fields given in (a 
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Schematic Model Configuration 



Fig. 5. Schematic of the regional QG model used in the assimilation 
study. The model streanifunciion M 7 and vorticily V 2 T are specified 
on the boundaries from climatological hydrographic data; during out- 
flow conditions the values for V : 4 7 are extrapolated onto the boundary 
from the interior model domain. Model mean quatities such as total 
kinetic energy were diagnosed in the interior dashed box with geo- 
graphical limits of about 29.5°-40.5°N. 28°-12°W. 

waves at the artificial boundaries. The damping is im- 
plemented by nudging the model fields toward the cli- 
matological background with a relaxation timescale lin- 
early increasing from infinity to 100 days over the outer 
10 grid points. The Geosat data are assimilated in the 
area enclosed by this 100-km wide friction layer. To 
avoid the influence of remaining boundary effects, in- 
tegral properties like kinetic energy or model-data misfit 
were evaluated in the inner dashed box, with geograph- 
ical extend from about 29.5° to 40.5°N, 28° to 12°W. 

3. Model results 

Starting from the initial model state given in Fig. 4, 
about two years of Geosat data were assimilated sub- 
sequently into the QG model. The initial and final times 
correspond to 15 August 1987 and 3 July 1989, re- 
spectively, giving a total period of 688 days. 

In order to study the sensitivity of the results to the 
specific choice of the relaxation coefficient, R l} 1 was 
modified between 0.5 and 5 days. A summary of those 
experiments is given in Table 2. Figure 6a illustrates the 
model history during experiment N21 from the total 
kinetic energy, T k = \H k ] j( V'P*) 2 dxdy, in all three layers 
(k = 1,2, 3). Also included in the figure is the observed 
kinetic energy T ohs . Noticeable are the large, up to 30%, 
variations in the observed energy with periods of max- 
imum and minimum energy coinciding with autumn and 
spring, respectively. The energy level in the surface lay- 
er decreases rapidly by about one-third during the first 
50 days, but the model tracks the observed energy vari- 


Tablk 2. Assimilation experiments and their parameters. Also 
shown are basin averages of model kinetic energy for all layers, as 
well as the mean rms model -data difference and cross correlation. 
Kinetic energy T in ergs per square centimeter. 


Expt 

T, T p 

(days) (days) T ilhs 

7', 

7\ 

T 

rmsA£ 

(cm) 

r 

NVH 

0.5 

1 1.54 

1.34 

0.43 

1.16 

2.4 

0.80 

Nil 

1 

1 

1.16 

0.33 

0.82 

2.6 

0.74 

N2I 

2 

1 

0.92 

0.24 

0.54 

3.0 

0.63 

N5I 

5 

1 

0.57 

0.15 

0.28 

3.5 

0.47 

N52 

5 

2 

0.80 

0.21 

0.21 

3.2 

0.59 

SD 

X 

— 

0.46 

0. 1 3 

0.19 

4.7 

0.02 


ations thereafter. The initial rapid decrease in T can be 
explained by the biharmonic friction that damps the en- 
ergy on small scales present in both the Geosat analysis 
and the initial model field of the surface layer. Based 
on the biharmonic friction coefficient a = 8 X 10 9 m 4 
s the damping timescale for the near grid scale signals 
is of the order of 25 days (Holland 1978). In contrast 
to the top layer, where real observations are used as the 
initial field, no adjustment process can be found in the 
lower layers. Instead, the variations in T 2 and 7^ are 
relatively small throughout the assimilation period for 



Fig. 6. (a) Total kinetic energy T k = i// A J J (V4 7 :* dx k ) dy in the 
three layers k obtained from experiment N21 as the spatial average 
within the inner dashed rectangle shown in Fig. 5. Also included in 
the figure is the observed kinetic energy T„ bs from Geosat data, (b) 
Area-averaged kinetic energy T ] of the surface layer that results from 
varying the relaxation coefficient between R = 1 day and R ' = 
5 days. The curve labeled Obs shows the observed energy similar to 
(a), and SD stems from the unforced reference run. Labels refer to 
experiments and parameters listed in Table 2. 






48 


JOURNAL OF PHYSICAL OCEANOGRAPHY 


Volume 27 


this choice of parameters. Fluctuations on the nudging 
timescale are subordinate, and a direct relation of vari- 
ations in T 2 and with those observed cannot be found. 

To illustrate the sensitivity of the model energy to the 
magnitude of the nudging parameter, Fig. 6b shows the 
surface layer kinetic energy T, from experiments with 
Rq 1 varying from 1 day to Note that the run, labeled 
“SD” (spin down), with R () 1 = ^ is driven only by the 
initial and lateral boundary conditions, with no data 
assimilated, and is provided here as an unforced ref- 
erence. 

As can be expected from Fig. 6b, the model energy 
decreases significantly with increasing relaxation time 
scale and for experiment N51 an energy level is ap- 
proached, which is equivalent to the unforced run SD. 
During this experiment the forcing of the nudging term 
is roughly balanced by its dissipative effect on basin 
average. The frictional character of the nudging term is 
illustrated in the figure from the initial decrease in T, 
which is significantly faster for all nudging experiments 
than for SD. Note the lag of the model behind observed 
variations in energy, which increases with decreasing 
coefficient R () . 

Figure 7 shows mean fields of (Figs. 7a-c) and FI 
(Figs. 7g-i) together with the variances of II (Figs. 7d- 
0, as they result from experiment N21 in all model 
layers, respectively (from top to bottom). The mean 
model state is in general agreement with the climatology 
but leads to smoother fields. Note that the near-surface 
frontal structure associated with the AC has become 
sharper through the data assimilation, and that in con- 
trast to the a priori assumption of a resting bottom layer, 
Fig. 7c shows a w estward Azores Undercurrent of about 
0.5 cm s '. The influence of the bottom topography on 
the abyssal mean circulation is especially clear in the 
western model domain where the W contours basically 
follow the bathymetry of the MAR. 

Maximum kinetic energy can be found along the AC 
in the top layer (Fig. 7d), with values continually de- 
creasing in downstream direction. The enhanced vari- 
ability in the upper left corner is related to an extension 
of the NAC, which is branching at that location into the 
Portugal Current. Contrary to the top layer, an increas- 
ingly homogeneous energy distribution in the meridi- 
onal direction is found in the subsurface. 

An analysis of the relative contribution of the various 
terms in the governing equations shows that the nudging 
term is dominating the model physics throughout the 
entire assimilation period with no apparent decreasing 
tendency. This implies an incompatibility between mod- 
el and data. The related discrepancy in eddy character- 
istics between the model and observations is apparent 
from the left column of Fig. 8 showing instantaneous 
fields of the streamfunction from N21 of all three model 
layers at day 210 after a one-half-year period of data 
assimilation (12 March 1988). Similar fields from the 
unforced run SD are shown in the right column of the 
figure for reference. With data assimilation, the upper 


layer is dominated by strong eddy activity with maxi- 
mum scales and amplitudes associated with the mean- 
dering AC frontal structure (Fig. 8a). A comparison w ith 
Fig. 4a indicates that eddy scales are significantly longer 
in the model than observed. Going from top to bottom, 
eddy amplitudes decrease and their scales increase. But, 
compared to the initial fields (Fig. 4), the eddy scales 
in the deep layers have been reduced through the con- 
tinuous data assimilation, and the pronounced initial 
vertical coherence has vanished. Without any assimi- 
lation, model fields are uncorrelated with those from the 
assimilation run soon after the initialization. Dominant 
scales are increasing with time and ultimately reach val- 
ues three times larger than those observed (see Fig. 9). 
In addition, the eddy amplitudes are decreasing signif- 
icantly, especially in the eastern part of the model do- 
main. 

A phase diagram (Hovmoller diagram) of the eddy 
field, which is given in Fig. 9 from 'P anomalies of both 
runs N21 and SD (relative to their 2-year mean) along 
35°N, indicates a predominantly westward propagation 
of eddy structures. In good agreement with the under- 
lying altimeter observations, phase velocities inferred 
from N21 are roughly 1.85 cm s“' in the surface layer 
(Fig. 9a). However, values are increasing toward the 
bottom. In the second layer (Fig. 9b), values of 2 cm 
s _l can be deduced. The westward propagating baro- 
clinic structures are broken up by much faster pertur- 
bations, which make it difficult to estimate phase ve- 
locities in the bottom layer (Fig. 9c). Those vertically 
coherent and fast fluctuations, which do not show up 
with continuous nudging on the full model domain (not 
shown), are dominant during several periods of low data 
coverage (compare Fig. 1 lc). Note that similar but much 
more vigorous fluctuations are apparent in the unforced 
case SD during the initial adjustment period of the mod- 
el after which the model has drifted into a completely 
different frequency-wavenumber regime than imposed 
by the observations. 

For a comparison of the model results with Geosat 
observations, Fig. 10 shows the observed sea surface 
elevation h obs (Fig. 10a) and its corresponding eddy sig- 
nal W (Fig. 10b), together with the resulting model fields 
from experiment N21 (Figs. 10c and lOd). Again, the 
model eddy fields are taken to be the anomaly relative 
to their 2-year mean and both fields correspond to 12 
March 1988 (day 210). Areas of large data uncertainty 
(with an estimated error exceeding 60% of the process 
variance) are hatched in Fig. 10a. They are found pre- 
dominantly in the eastern basin, where there is sub- 
stantial data loss. High amplitude anomalies present in 
both observations and the model fields are much dif- 
ferent in their spatial characteristics: while the obser- 
vations are dominated by anomalies with large merid- 
ional extent, the model tends to form isolated eddies. 
This difference can be understood from the fact that 
vorticity (f = V 2 ty) is assimilated rather than V by 
which larger scales are suppressed in the model results. 
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The root-mean-square difference between both fields is 
about 2.5 cm and the cross-correlation is r ~ 0.7. 

The rms difference and the cross-correlation between 
model and data eddy fields are given in Fig. 1 1 for the 
full assimilation period as a function of time. To indicate 
the range obtained from varying the relaxation timescale 
R 0 results from experiments Nil and N51 are dis- 
played. In the two cases, rms differences vary between 
2 and 4 cm and are significantly lower than an upper 
threshold limit of 6cm (labeled “Obs"), which indicates 
the differences of the observations relative to their initial 
field (persistence). Those latter values are reached asym- 
totically after the first 50 days of experiment SD, in- 
dicating the loss of correlation of the unconstrained run 
with the observations over that timescale. 

In order to give a rough estimate of the data coverage 
during the assimilation, Fig. 1 ic shows the number of 
valid daily Geosat observations over the model area. It 
is clear that periods of increased model errors and re- 
duced cross-correlation coincide with periods of sig- 
nificantly reduced data coverage. During those periods 
when up to 50% of the model grid is unforced by the 
altimetry, the fast transients, visible in Fig. 9, signifi- 
cantly degrade the agreements. A parameter T r = 1 day 
is sufficient to keep the cross-correlation as high as r 
= 0.8, but the increased value of i r — 5 days does not 
suffice to keep the model on track with the data during 
periods of massive data loss. 

4. Model-data comparison 

We proceed now with a test of the model simulations 
from experiment N21 against simultaneous and inde- 
pendent subsurface information available from hydro- 
graphic measurements and current meter data. Datasets 
used during the comparison are listed in Table 3. 

a. Hydrographic data 

Hydrographic data were measured on an eddy re- 
solving grid during two cruises in March 1988 (Kase et 
al. 1989) and in May-June 1989 (Hinrichsen et al. 
1993). The dynamic topography obtained from those 
data is used here as the observational counterpart of the 
model streamfunction. The data measured during March 

1988 were obtained during a 3 week period between 4 
and 26 March. Figures 12e,f show objective analyses 
of the resulting dynamic topography at 50 and 1000 
dbar relative to 3000 dbar, respectively (see Kase et al. 

1989 for details on the analysis). To allow a comparison 
with the model eddy field, the RBS climatology was 
subtracted from the hydrography. 

Shown in Figs. 12a,b is the eddy field from the Geosat 
analysis and its error map, taken from 12 March 1988 
(day 210). Figures 12c and 1 2d show the corresponding 
model eddy fields from layer 1 and layer 2, respectively. 
Discrepancies between the model simulations and the 
corresponding in situ measurements are, in particular. 


obvious close to the areas of enhanced observational 
errors. But otherwise the positions of various cyclonic 
and anticyclonic model eddies (labeled by capital and 
lowercase letters, respectively) and their scales are com- 
patible with the hydrography in both layers. It has been 
shown that the anticyclones A to D are the surface signal 
of subsurface meddies (Kase et al 1989; Stammer et al. 
1991). In contrast to the observations, however, feature 
B is separated from C by a trough and a second positive 
anomaly C' is connected to A and C by a ridge. Sim- 
ilarly, positions and scales of the simulated cyclones a- 
d agree with the observations. Note that in contrast to 
the Geosat analysis (Fig. 12a) the model is able to sim- 
ulate the observed positive anomaly B. In the second 
layer the amplitudes of the model eddies are signifi- 
cantly lower than have been observed, but their posi- 
tions agree. 

The correlation between the model fields and the hy- 
drography of r=0.58 and r=0.51 in both layers, re- 
spectively, is found to be rather low, but still significant 
on the 95% confidence level. (Our comparison is based 
on the objective analysis of the field data on a 50 km 
X 50 km grid. In order to evaluate the confidence test, 
we used only those 120 grid points that are close to the 
actual station positions. Due to the correlation scale of 
A = 100 km used during the hydrographic analysis, we 
considered only one-fourth of the grid points to be in- 
dependent, which leads to an effective degree of free- 
dom of 30.) Without data assimilation, no significant 
agreement of the model simulations with the observed 
hydrography is found (r < .1). 

For a description of the ocean state during the time 
of the second cruise see Hinrichsen et al. (1993). During 
that period, Geosat has lost substantial amounts of data 
in the model domain (compare Fig. 13b), and a direct 
comparison of both datasets is not possible. Now r the 
model dynamics have to carry all information into the 
data-sparse areas, which renders the comparison of the 
model results with the hydrography much more strin- 
gent. Results are shown in Fig. 13 with the model and 
altimeter fields taken from 21 May 1989. The hydrog- 
raphy resembles the dynamic topography 50 and 800 
dbar relative to 2000 dbar. As in Fig. 12, only the eddy 
signal is shown. Although not directly forced close to 
the hydrographic station net, the model fields generally 
resemble a spatially coherent picture of structures, 
which to some degree represent the field data. As an 
example, the features B and B' are simulated by the 
model in both layers, though with much reduced am- 
plitudes. The corresponding correlations between model 
fields and the hydrography is r=0.51 and r — 0.40, for 
the first and second layer, respectively. This is slightly 
below that during the 1988 survey, but still significant. 

h. Current meter data 

Time series of moored current meter data are available 
from the central model domain (see Fig. la) and allow 
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a comparison with velocity fields from all three model 
layers. Zenk et al. (1989) and Muller and Siedler (1992) 
provide a detailed description of the current meter data. 
Previously, these data have been used in a direct com- 


parison with Geosat geostrophic surface velocity esti- 
mates (Stammer et al. 1991; Stammer 1992). They will 
be used here to test the velocity in the subsurface model 
layers. Since the data from mooring “MW” were re- 
covered only above 1000-m depth, the following com- 
parison is based predominantly on data from mooring 
“Kiel276.” But some results from mooring “MW” are 
included in Table 4. 

Figure 14 shows the model velocities from experi- 
ment N21, subsampled once per day at the position of 
Kiel276 together with the observations. While the cur- 
rent meters at depth 1050 and 3000 m coincide with the 
mean depth of layer 2 and layer 3, the observations at 
450-m depth stem from a depth range that is rather close 
to the interface between the top two layers. Statistics of 
all time series are listed in Table 4. The mean model 
velocities are of the same order as the observed values, 
but the model kinetic energy is lower by a factor of 2- 
3 than the observed — much in agreement with previous 
results (compare Fig. 6). 

Despite the initial discrepancies between model and 
mooring velocities, which are likewise present between 
the mooring data and Geosat observations (Stammer 
1992), a reasonably good agreement is found between 
both datasets in the top two model layers, with observed 
velocity fluctuations being simulated by the model even 
on small time scales. A good example can be found in 
Figs. 14a,b around 10 June 1988. However, the model- 
data agreement is weak in the bottom layer. A similar 
result was found by Capotondi et al. (1995b) from their 
assimilation experiment in the Gulf Stream extension 
area. It was shown there that the model-simulated eddy 
signal is significantly correlated with current meter data 
down to about 1500-m depth but the agreement became 
insignificant below. 

Variance spectra of the velocity time series allow an 
analysis of frequency bands of enhanced agreement or 
disagreement. Shown in Fig. 15 are corresponding spec- 
tra of the v (northward) component. In layer 1 spectra 
from observed and simulated velocity time series agree 
on the 95% confidence level, and both follow roughly 
a —2 to —3 power law. Discrepancies between both are 
apparent only at frequencies above the nudging fre- 
quency (2 X 10 1 day '). On the contrary, the spectra 
from the deeper layers indicate systematic differences 
between the model and observations on all frequencies. 

In order to determine those parts of the spectra that 
are strongly affected by the surface constraints, results 
from the unforced experiment SD are included in the 
figure. (Similar to experiment N21, spectra from SD 
characterize only that data period covered by the moor- 
ing data; thus the initial adjustment processes are ex- 
cluded.) It follows that in all layers the model kinetic 
energy is significantly modified by the surface boundary 
conditions. In the surface layer a surplus of model en- 
ergy through assimilation reaches one to two orders of 
magnitude at all frequencies, while in layer 2 and layer 
3 the increase is found to be limited to frequencies below' 
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Fie.. 9. Phase diagram showing the model eddy field referenced to the 2-year model mean given in Figs. 7a to 7c and plotted 
along v = 1000 km (about 35°N) as a function of lime for all three model layers (top to bottom). Results are taken from 
experiment N2I in the left column (a to c) and positive areas are hatched. Panels (d) (ft show corresponding fields from 
experiment SD. 
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Fig. 11. (a) Rms model-data misfit as a function of time from 
experiments Nil and N51. Also included is the rms misfit of the 
observations with its initial field as a function of lime (labeled Obs). 
In (b) the cross-correlations from the corresponding fields are given. 
The curve labeled Obs in (b) shows the data persistence, that is, 
correlation of the observations with the initial field as a function of 
time, (c) Number of valid Geosat 1-Hz measurements in the model 
domain as a function of time. Note the pronounced correlation be- 
tween periods of increased misfit and reduced data coverage. 

4 X 10 2 cpd. In addition, a pronounced energy peak 
can be found near the nudging frequency in all layers. 
It is possible that the limitation of the energy increase 
to periods exceeding about 25 days is related to the 
dynamical response of the model to the bandpassed sur- 
face observations. It is likewise possible that this finding 
is associated with the baroclinic response of the model 
to the surface observations (Frankignoul and Muller 
1979). 

Although a visual comparison does not indicate an 
agreement between the two velocity time series in the 
bottom layer, a coherence analysis reveals a significant 


Table 3. List of available in situ datasets. 



Comment 


Hydrographic data 

4.3.-26.3.1988 

67 hydrographic stations on an eddy- 
resolving grid mean station spacing 
< 60 km (Kase et al. 1989) 

24.5.-20.6.1989 

Hydrographic measurements on an 
eddv-resolving station grid (Hin- 
richsen et al. 1993) 

Mooring data 

4.1 1.1987-8.1.1989 

Mooring MW (36°8.4'N. 18°23.4'W): 
lowpass filtered daily mean current 
measurements in 630 and 1034-m 
depth 

10. 1 1. 1987-10. 1 . 1989 

Mooring Kiel 276 (33°8.5'N\ 
21°57.6'W): lowpass filtered daily 
mean current measurements in 450/ 
630- m, 1050/1250-m and 3000-m 
depth (Zenk et al. 1989) 


coherence between model and mooring data in all layers 
at periods exceeding 80 days (Fig. 16). However, the 
analysis indicates also a peculiar phase lag, with the 
model leading the observations by 60° and 90° in layer 
2 and layer 3, respectively. This finding is puzzling but 
may indicate a systematic discrepancy between model 
and data: due to the low vertical resolution the baro- 
tropic mode presumably is getting too much energy dur- 
ing the data incorporation, which results in unrealisti- 
cally high phase speeds of induced eddy structures with 
relatively large scale. It is noteworthy that those phase 
lags are very sensitive to the presence of bottom to- 
pography. Without prescribed topographic structures, 
basically no coherence was found in the deep layers and 
topography with increasing amplitude leads to a de- 
creasing phase lag of the deep model relative to obser- 
vations. 

5. Summary and concluding remarks 

Perhaps the most important result from this study is 
that we find a significant correlation between simulated 
and observed structures of the ocean flow field in the 
eastern North Atlantic after constraining the model by 
Geosat surface SSH observations. This result is lending 
strong support to the hypothesis that even under realistic 
conditions at least aspects of the upper thermocline flow 
field of the ocean may be determined by combining sea 
surface observations with the physics of a numerical 
circulation model. Although studied extensively during 
various previous process studies using synthetic obser- 
vations, it was never possible to test the earlier en- 
couraging findings based on real data. The basic con- 
clusion that can be drawn from this study is that the 
results from process studies are indeed consistent with 
those found under realistic conditions. 

The assimilation of total surface elevation obtained 
from the sum of a climatological mean dynamic topog- 
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Fig. 12. (a) Geosat anomalies and model eddy field from the first (c) and second (d) layer, respectively, as they result 
for day 210 (12 March 1988). The hydrographic dynamic topography 50/3000 and 1000/3000 dbar is shown in panel 
(e) and (f), respectively, after the subtraction of the RBS climalologic dynamic topography fields. Positive values are 
hatched and cyclonic and anticyclonic eddy structures are labeled accordingly. In (b) the Geosat estimation error is given. 
Hatched areas indicate an estimation error exceeding 60% of the process variance. The cruise track is given in (e). 
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Fig. 13. As in Fig. 12 hut for the hydrography obtained during the period 21 May it) 20 June 1989. Shown is the dynamic 
topography (e) 20/2000 dbar and tf) 800/2000 dbar. respectively. Geosat and model results stem from 21 March 1989. 


raphy and Geosat SSH anomalies leads to a mean model 
state consistent with the initial a priori climatology in 
all model layers. It is shown that the simulated current 
field is consistent with hydrographic observations of the 
ocean state in the top two model layers. A significant 


coherence is also obtained between model velocities and 
time series from moored current meter data over most 
of the water column on time scales exceeding 80 days. 
The comparison with an unforced reference run indi- 
cates that the model is constrained by the surface ob- 
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Table 4. Comparison of model results with mooring data. Listed 
are mean values and standard deviations of the zona! (u) and merid- 
ional (v) velocity components from the moorings and corresponding 
model simulations. 




Mooring 



Model 
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a* 
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erf 
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tr 
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of 

Layer 

(cm) 

(cm 2 ) 

(cm) 

(cm 2 ) 

(cm) 

(cm 2 ) 

(cm) 

(cm 2 ) 




1 ) Mooring 1 

V1W 




SI 

1.3 

16 

0.5 

6 

3.6 

12 

0.9 

13 

S2 

0.6 

7 

0.4 

3 

0.6 

2 

0.2 

1 




2) Mooring Kiel 276 




SI 

0.4 

11 

-0.3 

16 

1.3 

8 

-1.3 

22 

S2 

-0.2 

3 

-0.3 

7 

0.4 

1 

0.1 

2 

S3 

-0.6 

1 

-0.2 

1 

-0.3 

1 

- 0.1 
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servations in all layers, leading to enhanced energy on 
small time and space scales. No agreement can be found 
between simulations and in situ data without the data 
constraint in both the surface and subsurface layers. 

In a parallel effort, Capotondi (1995a,b) assimilated 
Geosat altimetry into a QG model of the Gulf Stream 
extension area and compared resulting fields with var- 
ious types of simultaneous in situ observations. Al- 
though their focus was on the improvement of the cli- 
matological behavior of the mode! as a consequence of 
the assimilated mean surface observations, their con- 
clusions are in agreement with those presented here. In 
their study it was shown that the model-simulated eddy 
signal is significantly correlated with current meter data 
down to about 1500-m depth, but hardly below. It should 
be recalled, however, that due to the long adjustment 
time scale of the deep ocean potential vorticity, the con- 
vergence timescale toward the true state is much longer 
at depth than near the surface, even during twin exper- 
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Fig. 14. Vector time series of daily model velocities subsampled from experiment N21 at the position of mooring Kiel276. Also included 
are the corresponding vector time series from the mooring at 450 m, 1050 m, and 3000 m, respectively. 
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Fig. 16. (a) Coherence analysis between the iM northward) velocity 
component of the current meter data and model simulations given in 
Fig. 14. (b) Same as (a) but for the phase with observations leading 
for positive phase. In both panels thin solid, dash-dotted, and bold 
solid lines show r results from layer 1. 2. and 3. respectively. 


innents (e.g„ see HMR). Therefore, the poor agreement 
of the simulations here, and in Capotondi et al. (1995b), 
could partly be due to the relatively short duration of 
the experiments. 

Limits of present attempts in simulating the observed 
state of the ocean can be attributed to the poor quality 
of the Geosat altimeter data (observational errors and 
data distribution) and the limited representation of ocean 
dynamics in the regional QG model. Both shortcomings 
will be vastly improved in the near future. See, for ex- 
ample, Oschlies and Willebrand (1995) for a related 
attempt to assimilate Geosat data into the eddy-resolving 
Community Modeling Effort model of the North Atlan- 
tic. 

Presently, TOPEX/POSE1DON is providing sea sur- 
face height observations with an unprecedented error 
budget of less than 5 cm (Fu et al. 1994), and the pre- 
cision of those observations will even improve consid- 


Fig. 13. (a) to (c) Variance density spectra of the v (northward) 

velocity component given in Fig. 14 from the mooring data (thin solid (from top to bottom). Also included are similar curves resulting from 
line) and experiment N21 (bold solid) in all three layers, respectively the reference experiment SD (dashed line). 
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erably with increasing mission duration. Present-day 
ocean general circulation models, on the other hand, are 
run globally in realistic configuration on an eddy re- 
solving grid and differences with TOPEX/POSEIDON 
data are shown to be on the order of 10 cm on large 
scales (Stammer et al. 1996). Methods have been de- 
veloped that will allow the incorporation of altimeter 
observations into ocean general circulation models in 
an optimal yet practical way by accounting for errors 
in both data and models and by providing error estimates 
of the resulting fields (e.g. Stammer and Wunsch 1996). 
It can be anticipated that these methods will be used to 
obtain an estimate of the present ocean state by assim- 
ilating TOPEX/POSEIDON altimeter data jointly with 
information available from the diverse WOCE obser- 
vations. 
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